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Wo investigate the Hubbard model on the honeycomb lattice with intrinsic spin orbit interactions 
as a paradigm for two dimensional topological band insulators in the presence of interactions. Ap- 
plying a combination of Hartree Fock theory, slave rotor techniques, and topological arguments, 
we show that the topological band insulating phase persists up to quite strong interactions. Then 
we apply the slave-rotor mean-field theory and find a Mott transition at which the charge degrees 
of freedom become localized on the lattice sites. The spin degrees of freedom, however, are still 
described by the original Kane-Mele band structure. Gauge field effects in this region play an im- 
portant role. When the honeycomb layer is isolated then the spin sector becomes already unstable 
toward an easy plane Neel order. In contrast, if the honeycomb lattice is surrounded by extra 
"screening" layers with gapless spinous, then the system will support a fractionalized topological 
insulator phase with gapless spinons at the edges. For large interactions, we derive an effective spin 
Hamiltonian. 
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I. INTRODUCTION 

Topological insulators embody a new class of topo- 
logical states which have attracted great attention 
recently^^'* . The key for this flourishing development is 
the understanding that spin orbit interactions can real- 
ize topological insulating phases^"^. The theoretical pre- 
diction of such phases in real materials^°"^^ as well as 
their experimental observations^^"^^ are responsible for 
the success of this rapidly developing field. 

A topological insulator exhibits a bulk energy gap (like 
ordinary insulators) while the edge (or surface in three di- 
mensions) has gapless states which are protected by time 
reversal symmetry. The topological difference between 
a topological insulator and an ordinary band insulator 
is characterized by a Z2 invariant^ which is non-zero in 
the topological phase. The existence of this topologi- 
cal quantum number as well as the quantized spin Hall 
conductivity inspired the field^"^^^, in particular, it was 
shown that the topological insulator phase or in two di- 
mensions also called quantum, spin Hall (QSH) effect - is 
stable against weak disorder and weak interactions^°'^^. 
Inside a topological insualtor, Maxwell's laws of electro- 
magnetism are altered by an additional topological term 
with a quantized coefficient, which gives rise to interest- 
ing physical effects^^'^^. 

A major role has been played by a simple model intro- 
duced by Kane and Mele^'^ consisting of a hopping and 
an intrinsic spin orbit term on the honeycomb lattice. 
The Kane-Mele model (without the Rashba term) essen- 
tially consists of two copies with different sign for up and 
down spins of a model introduced earlier by Haldane^^. 
Haldane's pioneering work realizes the Quantum Hall 
effect without an external uniform magnetic field. It 
breaks, however, time reversal symmetry (necessary for 
the quantum Hall effect) which can be restored by tak- 
ing two copies with different signs for the spins together 
(as Kane and Mele did). Originally they proposed the 
model as realization of the QSH effect in graphene^'® 



and today it should be seen as a paradigm, a perfect 
theoretical model for topological insulator phases in two 
dimensions. The honeycomb lattice is definitely interest- 
ing on its own due to the striking development within the 
graphene community'^^ but also more exotic phenomena 
like zero modes^^'^^ have been discovered, for example. 
The honeycomb lattice has also attracted some attention 
in relation with exotic phases of light and the Jaynes- 
Cummings lattice model"^*. The additional spin orbit 
interactions now make the difference and are responsi- 
ble for the existence of a topological insulator phase^'^. 
Consequently, real materials with (strong) spin orbit in- 
teractions have been attracted particular notice^^""^^ . It 
was also shown that strong nearest- and next-nearest 
neighbor repulsions can imitate the intrinsic spin orbit 
interactions such that QSH phases are stabilized in the 
absence of spin orbit coupling*^"^^. Beside mercury tel- 
luride quantum wells and the Kane-Mele model, topolog- 
ical insulating phases in two dimensions are found to ex- 
ist in the Kagome lattice^^ and the decorated honeycomb 
lattice'*'' provided the presence of spin orbit interactions. 

Other aspects of topological insulators are disorder in- 
duced topological phases as predicted for the HgTe quan- 
tum wells''^''*'' and for three-dimensional systems^'^'''''^ 
and the proposed existence of axions on the surface of 
bismuth-tin alloys^^. Axions were postulated more than 
30 years ago in the context of the standard model^^ and 
their effective action has now been recovered in topologi- 
cal insulators raising hope to detect this dynamical axion 
field experimentally. Moreover, a QSH phase in ferro- 
magnetic graphene was predicted^^ which is protected 
by the product of charge-conjugation and time reversal 
symmetry. Most recently, a new family of topological in- 
sulators has been discovered^*'^^ in ternary Heusler com- 
pounds. Their additional open /-shell element might be 
the key for the realization of exotic topological effects. 

Another promising path for the realization of topo- 
logical phases and, in particular, QSH phases consists 
of cold atomic gases loaded into optical lattices^^ which 



are subjected by a synthetic magnetic field. Such a field 
has a similar effect on the neutral atoms as a magnetic 
field coupled to electrons and has been demonstrated 
experimentally^^. A possible experiment to realize a 
topological insulator was proposed recently^^"^^. In this 
spirit, a realization of a topological band insulator seems 
to be feasible in the near future with possibly two major 
advantages: (i) tuning of the topological insulator band 
gap or of the details of the engineered Hamiltonian and 
(ii) availability of onsite-interactions (Hubbard model) 
with tunable interaction strength. 

In this paper, we investigate the Hubbard model with 
intrinsic spin orbit interactions on the honeycomb lat- 
tice which corresponds to the Kane-Mele (KM) model 
with interactions. Some aspects of the interacting KM 
model was studied in Refs. 62,63. A general theory of 
interaction effects in topological insulators has been pro- 
posed introducing a topological order parameter in terms 
of the full Green's function^''. We consider the half-filled 
case at zero temperature. While the (non-interacting) 
KM model is known to realize a topological band insula- 
tor (TBI) phase, it is also expected that for sufficiently 
strong electron-electron interactions magnetic order will 
take place. Therefore we want to clarify what happens 
and which phases are present when adding interactions 
- ranging from very weak to very strong. We focus on 
the dominant phases at finite spin orbit coupling. For 
very weak or no spin orbit coupling additional (spin liq- 
uid) phases might exist but are beyond the scope of this 
paper. First, we show that interestingly the TBI phase 
subsists up to quite strong interactions. 

Then, applying the slave rotor mean-field procedure, 
we investigate the limit of stronger interactions where the 
charge degrees of freedom form a Mott insulator whereas 
the spin degrees of freedom are described by a renormal- 
ized KM model. In a Mott phase, adding a particle at 
a given site costs the Hubbard onsite energy U (in con- 
trast, excitations carry a well-defined momentum in the 
TBI phase). At the mean-field level, this phase has all the 
properties of a spin liquid (which preserves time-reversal 
symmetry) with gapless spinon excitations at the edges 
and is characterized by a hidden order parameter in the 
spin sector similar to that in the original KM model^'^'^. 
On the other hand, one should not underestimate the ef- 
fect of dynamical compact U(l) gauge fields, especially 
in two dimensions'^. 

In particular, one predicts^' that such a spin liquid 
phase (with gapless edge spinous) found at the mean- 
field level can only be stable beyond the mean-field limit 
if other gapless layers (spinous) are present to screen the 
gauge field and suppress the gauge fluctuations. Poten- 
tial candidates can be found in Refs. 67,68. Furthermore, 
Mott physics will suppress the single-particle tunneling 
at the edges'^ such that the lowest relevant coupling 
between layers is the usual spin-spin interaction which 
may remain irrelevant®', then preserving the gapless edge 
spinous. Phases exhibiting similar spin-charge separation 
were also reported in other systems^"''^^''' and for topo- 
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FIG. 1; (Color online) Phase diagram of the isolated hon- 
eycomb layer where the proliferation of instantons produces 
a Neel order in the XY plane already in the entrance of the 
Mott phase. Above the red dashed line the SDW phase can 
be described in terms of a mean-field Hartree-Fock theory 
whereas below the red dashed line the easy plane Neel order 
emerges as a result of subtle gauge fluctuations beyond the 
mean-field solution. (The precise nature of the "transition" 
associated with the red dashed line is beyond the scope of this 
paper.) 
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FIG. 2: (Color online) Phase diagram in the presence of ad- 
ditional screening layers (with gapless spinous) allowing to 
screen the gauge field and therefore stabilize the Fraction- 
alized TI phase found at the mean-field level. Here, charge 
degrees of freedom are in the Mott regime and spin degrees 
of freedom form a spin liquid with gapless edge spinous. 



logical insulators in the presence of a tt fiux^'^'^^. In con- 
trast, if the honeycomb layer is isolated then the prolif- 
eration of instantons will fatally result in a Neel ordering 
in the XY plane^^'^'^. The two distinct scenarios at the 
Mott transition are reported in Fig. 1 and 2. Following 
Ref. 66, Fractionalized TI refers to the spin liquid type 
Mott phase with gapless spinous which preserves time- 
reversal symmetry and SDW in the two figures always 
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refers to the occurrence of a spin density wave formed in 
the XY plane. 

For very large interactions, applying a conventional 
Hartree-Fock procedure and deriving an effective spin 
Hamiltonian, we show that SDW phases with XY order- 
ing are allowed on the honeycomb lattice when adding 
the spin-orbit term coupling next nearest neighbors. 

The paper is organized as follows. In Sec. II we intro- 
duce the KM model, re-derive some of its basic proper- 
ties, and introduce the (Hubbard) interaction we consider 
throughout the paper. In Sec. Ill we apply the Hartree 
Fock method in order to show that a conventional SDW 
phase with ordering in the XY plane appears at large 
U. In addition, we derive an effective spin model. Then, 
in Sec. IV, we use a mean field approach in momentum 
space as well as the slave rotor picture to argue that the 
TBI phase as present in the original KM model is stable 
beyond renormalization group results®'^^ up to moderate 
interactions. Then, we apply the slave rotor theory of 
Florens and Georges^^"^'' and discuss the intermediate 
region and the gauge field effects more thoroughly. 



II. MODEL AND GENERAL 
CONSIDERATIONS 

The Kane-Mele (KM) model'^^^ which might be con- 
sidered as a spinful version of the Haldane model consists 
of two parts, a nearest neighbor hopping term and a sec- 
ond neighbor hopping spin orbit term on the honeycomb 
lattice, 
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FIG. 3: (Color online) Left: Honeycomb lattice consisting of 
two interpenetrating triangular lattices, A (dark blue dots) 
and B (cyan dots), with its lattice vectors ai and 02 (dashed 
arrows). In addition, the three nearest neighbor vectors 6i 
{i = 1,2,3) are shown connecting the two sublattices (solid 
arrows). Right: Corresponding Brillouin zone with the two 
inequivalent Dirac cones K and K' and the high-symmetry 
points r and M. 



and ^5 6== ±(02 ^0,1)- In what follows we set a = h = 1. 
Throughout the paper N\ denotes the number of unit 
cells, while N is the number of particles. Hence, the num- 
ber of lattice sites is 2A^a and at half filling TV = 2A'a. 
If needed, we refer to the sublattices A and B as Aa and 
Ab- The previous definitions imply 
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As a first step we wish to reproduce the energy bands 
due to nearest neighbor hopping and switch to momen- 
tum space via 



Here Cj^ is an electron annihilation operator either on 
sublattice A or B (then denoted by or bi„, respec- 
tively) fulfilling the fermionic standard anti-commutation 
relations {cjg.,c|^,} — 6ij6cra'- As usual t is the hopping 
integral and A is the spin orbit coupling, (ij) denotes 
nearest neighbor and <^ ij 3> next nearest neighbor sites, 
is the third Pauli matrix and i^ij = ±1 as discussed 
below. Throughout the paper we consider the Rashba 
spin orbit interaction to be zero. The lattice vectors of 
the honeycomb lattice are given by 
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and shown in Fig. 3. The lattice vectors have the length 
-\/3a while the lattice spacing a is the distance between 
neighboring atoms A and B. Note that our notation of 
the honeycomb lattice is adapted from the review of Cas- 
tro Neto et al.^^. We further have the nearest neighbor 
vectors 
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^(1,-V3), 



<53 = a(-l,0) (3) 



which are also shown in Fig. 3. The six next-nearest 
neighbor vectors S[ are given hy 6[ 2 — io^i, '^34 = ifl2, 
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The function g is given by g = g{k) ~ t X]^=i e''^^^ . Here 
we used the unitary transformation matrix 
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FIG. 4: (Color online) Left: Flux configuration per spin and 
sublattice associated with the intrinsic spin orbit term. Right: 
The flux configuration of one sublattice corresponds to a stag- 
gered triangular flux lattice. 



to diagonalize Hk via T^IIkTo — diag(— Cal- 
culating ±|g| explicitly results in the well known tight 
binding spectrum of the honeycomb lattice 
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two particle hole symmetric bands which touch each 
other at the six corners of the Brillouin zone (BZ) corre- 
sponding to two inequivalent points. Expanding around 
these special points reveals a linear dispersion which gives 
rise to the name Dirac points. The positions of two in- 
equivalent Dirac points which are shown in Fig. 3 are 
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Although it is very convenient to expand around the 
Dirac points and formulate a Dirac equation on the hon- 
eycomb lattice we will keep throughout the paper the full 
tight-binding model. 

As a second step, we consider the intrinsic spin orbit 
term*^ of the KM Hamiltonian (1). The expression Vij 
gives ±1 depending on the orientation of the sites. A 
formal definition is 
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where di and (^2 are the unit vectors connecting the sites 
j and i. Essentially, making a left turn yields "—1" while 
a right turn "-|-1". Note that hopping from a site of sub- 
lattice A in direction Sj would yield the opposite sign 
than hopping from a site of sublattice B in the same di- 
rection. Hence we should keep in mind that Vij oc 
(where is again the third Pauli matrix). As we will 
see below the spin orbit term opens a gap in the bulk. For 
completeness, we notice that other possible terms which 
open a gap in the spectrum are different from the spin 
orbit term. Such other terms, like a staggered sublattice 
potential Hst = X)i ^i'^la'^ia where = 1 on sublattice A 



and = — 1 on sublattice B, result in an ordinary band 
insulator and not in a topological phase since the gap 
is spin-independent. The spin orbit term preserves the 
original unit cell. We have shown the corresponding flux 
configuration per spin and per sublattice in Fig. 4. The 
net magnetic flux through a plaquette is zero following 
Haldane's idea^"'. The flux pattern for one of the sublat- 
tices corresponds to a triangular staggered flux lattice. 
Transforming the spin orbit term to momentum space 
leads to: 
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where ^'^ = {a\^^h\^,a}f^^,b\^^ and cr^r^ is understood 
as a 4 X 4 matrix, ct^t^ — diag(l, —1, —1, 1) (cr^ for spin 
and for sublattices). 

The bands of the KM model are now obtained by 
diagonalizing the 4x4 matrix of "H — Ht + V-so — 
^k^^k'^k with the matrix 
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Blank entries should be thought as zeros. As the matrix 
consists of two independent 2x2 matrices, the diagonal- 
ization procedure is identical to Eq. (8) when replacing 
Tg by T-^ and T^. The exact form of the transformation 
matrices, and T^, is explicitly given in Sec. IV. The 
single particle spectrum of the KM model in the infinite 
system consists of two two-fold degenerate energy bands 
(reflecting the Kramers degeneracy), 
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which are plotted for several values of A in Fig. 5. The 
spectrum in Eq. (15) is still particle-hole symmetric. An 
important feature is that an infinitesimal value of A opens 
an infinitesimal gap at the Dirac points. For evaluating 
the gap size due to the spin orbit term we know that only 
the Dirac points K, K' as well as the zero-energy lines in 
7 play a role and it is, hence, sufficient to consider these 
special points. At the Dirac points, we find 
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FIG. 5: (Color online) Energy bands {t = 1) of the Kane- 
Mele model for (a) A = 0.05, (b) A = 0.2, (c) A = 0.5, and 
(d) A = 1.0. The "path" through the Brillouin zone is taken 
as shown in the inset of (a). 

At the zero-energy lines of 7, we find (without loss of 
generality we consider here the line ky = only) 

e(fc., 0) = ^t^ (5 + 4cos(3fc,/2)) > \t\ . (17) 

The minimal value t is reached at the M point of the 
Brillouin zone. We summarize that the dispersion grows 
linearly with A at the Dirac points, but immediately when 
the value of A = l/(3-\/3)i is reached, the gap remains 
constant with a gap size 

A = 2t (A > 0.193 i) . (18) 

We further consider an ordinary next-nearest neighbor 
hopping term without spin-orbit interaction for reasons 
which will become clear in Sec. V. This term is identical 
to 'Hso when omitting i, A, and Vij and replacing a^^, 
by Sa-a-'- Hence we will find a function 172 instead of 7, 

52(fc) = 2cos(\/3fcj,) + 4cos(\/3fcj,/2)cos(3fc:,/2) . (19) 



In contrast to the spin orbit term, it breaks particle-hole 
symmetry but does not open a gap at the Dirac points. 
As already mentioned this term is not present in the KM 
model (1) but will become relevant in Sec. V. 

The aim of this paper is to investigate the effect of a lo- 
cal Hubbard interaction to the KM model: the Hubbard 
term reads 

H'j^U^n^^n,^ . (20) 

i 

We consider the case of half filling which allows us in 
principal to rewrite the Hubbard interaction as 




This is a particularly convenient formulation for the slave 
rotor theory, we will, however, use the form (20) as well. 
Note that Eqs. (20) and (21) are identical at half filling. 



III. SDW PHASE AT LARGE U 
A. Mean field arguments 

In the past, a transition from the semi-metal (Dirac liq- 
uid) to a SDW phase has been evidenced in the context of 
the ordinary Hubbard model (A — 0) on the honeycomb 
lattice, when increasing the strength of the onsite inter- 
action. Within the Hartree-Fock approach, this transi- 
tion takes place^^ at Uc — 2.23 1. Within quantum Monte 
Carlo (QMC) the transition was found''® at Uc ~ 4.5 — 5 i, 
and by means of dynamical mean field theory (DMFT) 
the transition''^ occurs even for Uc > 10 1. While the 
critical value of the interaction parameter strongly de- 
pends on the used method there is no doubt about the 
existence of the SDW phase for strong interactions. The 
reason for the occurrence of a SDW phase is simply the 
bipartite nature of the honeycomb lattice. 

In this Section, we will apply the Hartree Fock method 
to determine the phase boundary Uc{X) at which it be- 
comes favorable to decouple the Hubbard interaction 
(20) in terms of the sublattice magnetizations mf — 

iEAaUAb 
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where rn = rii^ + 714, nf^ = al^a^^, and UJ2in'^/4' = c 
is a constant in the SDW phase. The SDW phase at 
large U imphes a Mott insulating phase for which we can 
evaluate (nf) explicitly. Indeed, for U — ?> +00, the Mott 
state is described by the exact wavefunction 



\M) = l[albl\0) , 



(23) 
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where a is either f or l while a "points" in the opposite 
direction. We find 



(Af I nf^ + n\ + 2^4714 |A/) = 1 . 



(24) 



Only one of the first two terms contributes depending on 
the sublattice site i belongs to; the third term is always 
zero due to the definition of \M). Hence (n|) = 1 and c 
is a constant in the SDW phase. We further assume that 
TTif^^^^ — rrS^/^K In momentum space the mean field 
(MF) decoupled Hubbard interaction reads 
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While we could keep and independently, here we 
will search only for an antiferromagnetic solution. To be 
more precise we are considering only an Ising-like order 
parameter. In principal, one could also treat the full spin- 
rotational problem {e.g. within a a model representation 
), we expect, however, no fundamental discrepancies with 
the simpler approach used here. In order to find the SDW 
phase we set 



(26) 



Hence the mean field Hamiltonian can be written as 
^HF ^ jj^i2 diag(-l, 1, 1, -l))*fc 



U 2 



(27) 



We notice that the mean field Hamiltonian coincides with 
the original KM model when replacing 7 by 7 — [/to/2 
up to additional constants. Now we write the mean field 
free energy at T = as 



= E ( - 2V IffI' + - Um/2f) + 



Um^NA 



(28) 

Minimizing the free-energy with respect to m yields the 
following self-consistent equation. 



mU/2 - 7 
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FIG. 6: (Color online) The numerical solution of Eq. (29) is 
shown {t — 1). For A = 0, we confirm tliat the SDW transition 
occurs at Uc = 2.23 1 in agreement with the result of Sorella 
and Tosatti^*. With increasing A, Uc increases up to 8.55 at 
A = 1.0. The fact that Uc{X) increases at finite A may be 
understood from the effective spin model. 



The solution of Eq. (29) provides the phase boundary 
as shown in Fig. 6. For A = we reproduce the earlier 
result by Sorella and TosattiJ* which is Uc = 2.23 1. With 
increasing A we find that Uc{X) increases. The reader 
may notice that formally we should add a Fermi function 
in Eqs. (28) and (29) since it may play a role for larger 
values of to. Here we focus on the phase boundary only 
which implies small values of to and neglect this point. 

For the usual Hubbard model (i.e., A = 0) on the hon- 
eycomb lattice the occurrence of the SDW phase is not 
surprising since we know that on any bipartite lattice 
an antiferromagnetic insulator, the SDW phase, will be 
favored in the limit of large U (at least at half filling). 
The effective Hamiltonian which describes the low-energy 
behavior of the Hubbard model for large values of U/t 
is the antiferromagnetic Heisenberg model. It is, how- 
ever, rather unclear whether the Mott transition, i.e., 
the phase transition from a semi metal into a gapped in- 
sulator phase, occurs simultaneously with the transition 
from a disordered spin state into the SDW phase. This 
unclearness is reflected in a current debate^""*^ . 

While in the early works^^ no indication was found for 
two separate phase transitions (Mott-Hubbard and mag- 
netic phase transition), Lee and Lee reported a possible 
realization of the nodal spin liquid^° directly at the Mott 
transition. While Herbut favored a direct semi metal- 
SDW transition using a large N approach^^, Hermele^^ 
proposed the stability of the SU(2) algebraic spin liquid 
in a small region followed by a valence bond solid phase. 
A recent QMC investigation^^ rather predicts the exis- 
tence of a resonating valence bond (RVB) phase between 
the Dirac-liquid and the SDW phase. 

Within the Hartree Fock procedure presented above we 
cannot answer the question, since the SDW order already 
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implies the Mott phase and docs not tell anything about 
the question where the phase transition into the Mott 
phase or in a spin liquid phase occurs. It is, however, 
clear that a possible spin liquid phase must be somewhere 
below the transition Uc = 2.23 1. The SDW phase is 
definitely the upper boundary of such a scenario. The 
fact that C/c(A) increases at finite A may be understood 
from the effective spin model which we will discuss now. 

B. Effective spin model 

In what follows we will investigate the behavior of the 
spin orbit term in the strong coupling limit U ^ oo 
with t and A keeping fixed. Similarly to the usual Hub- 
bard model, we expand the Hamiltonian in powers of 
t/U. The spin model can be derived in a systematic 
way as shown e.g. in Ref. 86, but essentially we have 
to consider the second order process of the Hamiltonian 
^so — S<ij> (^so)ij with the additional prefactor 
—2/U. The minus sign respects second order perturba- 
tion theory which always lowers the energy: 

('Hso)ij ('Hso)ji 
= -UijVjiX^ (al^a.-t - J ("jt'^iT ' sV^J (30) 

= A^ {2SfS] + 2S^Sy - 25f 5| - ^mnj + n,^ . 

Without loss of generality we have considered the hop- 
ping process on sublattice A, but there is no difference 
with the equivalent process on sublattice B. At half fill- 
ing and for U — >■ 00, we can assume rii = nii^ + riii ~ ^ 
and neglect the last terms which are constant. Together 
with the mentioned factor —2/U we find the efi'ective spin 
model 

a = I J2I {-S^S^ - SfS] + S^S^) (31) 

with the exchange coupling J2 = 4A^/C/ on a triangu- 
lar lattice (where we assume that the sum counts every 
nearest neighbor pair only once). The spin model then 
consists of a XY term which favors ferromagnetic order 
and a Z term which favors antiparallel alignment of the 
spins. From the nearest neighbor hopping term we ob- 
tain an isotropic antifcrromagnctic Hciscnbcrg term with 
Ji = At'^/U. This term stabilizes the antifcrromagnctic 
order, i.e., the SDW. The J2 term (more precisely, its 
^-component) competes with the Ji term. This tends to 
explain the increase of the critical interaction to reach 
the SDW order in the phase diagram Fig. 6. 

On the other hand, the a;r/-component of the J2 term 
favors ferromagnetic order in the XY plane on each sub- 
lattice. While the ordering vector might point in any di- 
rection when A = 0, we assume that ordering within the 
XY plane is preferred immediately when A ^ 0. This can 
also be seen from energetic arguments. Once the ground 



state is XY ordered, "f and "|" in |M) refer to the x- 
component of spin. In other words, acts like a spin 
flip operator. Consequently, any operator containing Sf 
has a zero expectation value; it particularly implies that 
for U +00, (M| Efe 7*ilf^^T^*fc \M) = 0. 

This suggests that the SDW phase (with preferable XY 
order) might even persist for J2 > Ji which is beyond the 
scope of this paper. We conclude this part by making a 
brief comparison with the ordinary J1-J2 model on the 
honeycomb lattice with J1/2 > 0. For weak values of 
J2/J1, a Neel order is present. For moderate values of 
J2/J1 ~ 0.3 a resonating valence bond (RVB) phase was 
proposed^^. For stronger frustrations, a valence bond 
crystal exhibiting a spin gap is reasonable^^. 



IV. STABILITY OF TBI PHASE 

In this Section, we provide several arguments estab- 
lishing that the topological band insulator phase present 
for A > in the original KM model is stable (against 
Mott physics) not only for weak intcractions^'^'^^'^^ but 
also for moderate interactions U ^ t. We first consider 
a mean field approach in momentum space which shows 
that the effect of ?7 < 2 1 for A > 0.2 1 does not affect 
the insulator phase and, hence, should be irrelevant for 
the topological band insulator phase. Then we intro- 
duce the slave rotor theory of Florens and Georges"^^'^^ 
and argue that this provides a rigorous proof concerning 
the stability of the TBI phase beyond the perturbative 
regime. We will discuss both approaches to show that 
the obtained results are valid beyond the renormaliza- 
tion group method^. In addition, wc will briefly discuss 
the Z2 topological invariant for the present situation. 



A. Mean field arguments 

The tight binding approach starts - after Fourier trans- 
formation with the momentum operators for both sub- 
lattices (ofeo-j bka)- Diagonalization of the tight binding 
matrix Hk introduces a new set of operators associated 
with the bands. As the spin remains a good quantum 
number (see e.g. Eq. (14)), we call the new operators Uka 
and Ika where u and I refer to upper and lower band. 

The explicit transformation matrices between the sub- 
lattice basis {uka^bka) and the band basis {lka,Uka) are 
given by 
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We define the functions 



a± = a±{k) = Af±d± , 



(34) 



term, as we will sec, where the chemical potential is given 
by t//2. The whole procedure does not use any additional 
assumptions and works for any interaction strength U. 



where 



and 



I3± = P±{k)=Af± , 



A4 



(35) 



(36) 



(37) 



Note that g, 7, e and then also d±, a±, j3±, and N± are 
fc-dependent. But for the sake of clarity we omit the k- 
dependence. To give the reader a better idea about a± 
and /3±, we have plotted |a_(fe)p for A/t — 1.0 and for 
X/t = 0.2 ill Fig. 7. Let us mention that d± is a complex 
function and hence a±, while (3± is real. Nonetheless, 
to make the equations more symmetric, we will use the 
complex conjugate of /3± as well. Technical aspects and 
mathematical considerations associated with the change 
of basis are presented in Appendix A. There, we show 
useful formulas like \a±\'^ = \l3zf\'^ and obtain the impor- 
tant result: 



^ |a±(fe)r= J2 l/3±(fe)l' = 



Na 



(38) 



fee BZ 



ke BZ 



Now we want to investigate the effect of the Hub- 
bard term on the topological band insulator state more 
deeply. For that purpose, we transform the interaction 
term into the band basis. Assuming A > 0.2 1 we know 
from Eq. (18) that the gap size A = 2 1 is large and hence 
we can neglect all terms containing operators of the up- 
per band (since it costs roughly the energy A to make 
any process between lower and upper band). Then we 
decompose the remaining term in a standard way. Even- 
tually the Hubbard term reduces to a chemical potential 



53 '^k+qt'^k'-ql^k'l'^kt + ^k+qt^k'-ql^k'l^kt 
kk'q 



J2 {Ak,k',q) + B{k,k',q)) 



I I 



kk'q 



k+q'[ k' — q^ k' ^ k\ 

(39) 



where we have suppressed all terms containing operators 
Ukcr or in the last line. The prefactors A and B are 
given by 



A{k, k',q) = a*_ (fe + q)al{k' - q)a+ {k')a- (fe) 



(40) 



(41) 



B{k,k',q) =0*_{k + q)Pl{k' - q)l3+{k')p.{k) . 
Now the Hubbard term will be decomposed as follows: 

^k+qt^k'-qi^k'l^kt ~ ^loi^kt^kt'^^k'i^k'l + 

+ ^kf^kt^Qoilk'llk'i) ~ ^<l0{lkt^kt)i^k'i^k'i) 

Then the prefactors reduce to 

A{k,k',q)Sqo=\a-{k)\^\a+{k')\^ , 
B{k,k',q)6qo=\l3.{k)f\0+{k')f . 



(42) 



We find the MF-decoupled interaction term in the new 
basis where we assume (nj^^) = = 1 for ct =t,-i 

since the lower band is completely filled. Therefore, we 
get: 



~ E ]^ (^■^(^' + Ak', k, 0)ni-f + B{k, k', 0)4^ + B{k', k, 0)4^ - B{k, fe', 0) - ^(fe, fe', 0) 



+ |^_(fe)P^^|/?+(fe')|Xt + l^+Wr]^El/5-(*=')lX4) - UNa/2 



(43) 



u 
'2" 



( (|a_(fe)|2 + |^_(fe)H ni^ + (|a+(fe)|2 + |/3+(fe)|2) 4j = Yl < ' 



with the effective chemical potential p.^^ = U/2 and n'^, = n^.^ -|- nj^j^. In the last line of Eq. (43), we have dropped 




FIG. 7: (Color online) Exemplarily the function jQf_(fc)p is 
shown in the BZ for \/t = 0.2 (left) and \/t = 1.0 (right). 
The dark blue region corresponds to a value of zero, while 
the white region to a value of 1. The function |a+(fc)p is 
identical to |Q_(fc)p by interchange of white and dark blue 
regions. The black lines mark the boundary of the BZ and 
the black dot in the center the F point. 

the constant term. For A > 0.2i, as the gap always spans 
from —t to i, we conclude that as long as U < 2t the 
effective chemical potential lies in the gap, i.e., the U 
term does not affect the TBI gap. Thus we can assume 
that the physics in the KM model will be unchanged 
and consequently the edge modes persist (they will be 
described by the helical Luttinger liquid theory as a result 
of interactions^''). 

B. General slave rotor arguments 

Now we will apply the slave rotor formalism which al- 
lows to address correlated electron systems at weak up to 
moderate interactions. The method has been introduced 
by Florens and Georges^^'^^; we also refer the reader to 
the review by Zhao and Paramekanti^^ . 

Within this approach^^"^^ , the original fermion opera- 
tors will be rewritten by a product of a fermionic operator 
the spinon or auxiliary fermion, and a phase factor 
e'^' , the rotor, 

c^a = e^'' U ■ (44) 

The idea is that the original fermions Cg. are represented 
by a collective phase degree of freedom 9 (conjugate to 
the total charge) and auxiliary fermions fa-. Introducing 
an additional variable, the angular momentum L cx idg 
associated with a quantum 0(2) rotor 9, simplifies then 
the original quartic interaction between the fermions as 
it is replaced by a simple kinetic term oc L^. State 
vectors in the new Hilbert space should have the form 
1 5*) = l^y) l^'e). The price we have to pay for the whole 
rewriting of the original problem is that the new Hilbert 
space is enlarged compared to the original one since un- 
physical states are present. To resolve this problem we 
have to impose a constraint, 

J2flf,^ + L, = l. (45) 



Since the original fermion operators fulfill anticommuta- 
tion relations also the spinon (or auxiliary fermion) op- 
erators do so. The reader may notice that in the rotor 
condensate phase the original electron and spinon op- 
erators are proportional, and one will find a situation 
where the spinon band structure describes physical elec- 
trons. In contrast, when the rotor is uncondensed, there 
is spin-charge separation and the spinous are emergent 
charge-neutral quasiparticles carrying spin only. The 
term "spinon" should not imply that the particles as- 
sociated with the new / operators obey fractional statis- 
tics in the spirit of the elementary excitation of the one- 
dimensional Heisenberg antiferromagnet. 
Rewriting the hopping term Ht yields 

nt = -^EE (/-^/?^ e-^'^ +h.c.) (46) 

where refers to sublattice A/B and 9ij = 9i — 9j. 
The spin orbit term "Hso has a similar form, 

nso = *A E E'^^^<-'/-/.<^'e"^''^ ■ (47) 
Rewriting the Hubbard term (21) which is local yields 

^^ = |EfE-f.-i)'-|E^?' (48) 

where we have used the constraint and = fl^fia- 

Without proceeding further, the introduced formalism 
already allows to read off the following results. The Hub- 
bard interaction U affects the rotor sector only. As long 
as U < t the rotors will condense favoring t he uniform 
ansatz 9ij = 0. It implies that the auxiliary fermions fia 
behave like the original electrons since exp {zti9ij ) = 1 
far away from the phase transition. Superfluid or Bose- 
condensed phases are known to be robust (roughly up 
to U ^ t). We can assume, hence, that the superfluid 
phase of the rotors also persists against moderate in- 
teractions before the phase transition in the insulating 
phase occurs. Both approaches presented in this section 
are beyond renormalization group (RG) arguments (the 
perturbative regime). Here, we have chosen to rewrite 
the Hubbard term in the rotor variables. 

In the seminal paper by Kane and Mele^ the stabil- 
ity of the TBI phase is briefly discussed. The derived 
RG equations indicate that additional Goulomb interac- 
tions increase the spin orbit gap size and does not destroy 
the TBI phase. The RG procedure, however, is only ap- 
plicable in the perturbative regime where t ^ U, X. In 
contrast, the arguments presented here provide evidence 
that for A > 0.2 1 the TBI phase is stable up to C/ ~ t 
(against the Mott phase) reaching the strongly interact- 
ing regime. As a last point we shall mention that the 
shown stability of the TBI phase implies also the stabil- 
ity of the edge modes (which are described by the helical 
Luttinger liquid theory as a result of interactions^^). 
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FIG. 8; (Color online) Brillouin zone with the four time- 
reversal invariant momenta F and M. In addition, the re- 
ciprocal lattice vectors (with half of their length) are shown 
(red arrows). 



In Appendix B we present the slave rotor approach us- 
ing a simple approximation to decompose Eq. (46). We 
will restrict ourselves to the case A = 0. When perform- 
ing the mean field procedure it turns out that the used 
approximation leads to results which are not so reliable 
in finite dimensions. Therefore, in Sec. V, we will apply 
the cr model representation of the slave rotor theory to 
find the transition for finite spin orbit coupling where the 
rotors undergo a quantum phase transition from super- 
fluid to Mott insulating phase. Before we address this 
issue, we briefly discuss the stability of the TBI phase in 
the context of Z2 invariants. 



C. Z2 invariants 

For the (non-interacting) KM model considered here, 
there exists in principle two non-trivial topological in- 
variants, the spin-Chern number of Sheng et al?^ and 
the Z2 invariant proposed by Kane and Mele^'^. 

The spin-Chern number can be evaluated by integrat- 
ing the Berry curvature of a fiber bundle obtained by 
imposing twisted boundary conditions'^. The proce- 
dure demonstrated in Ref. 21 implied that the spin Chern 
number is a robust topological invariant. Essentially the 
idea is that the system conserving decouples into two 
independent Hamiltonians for the up and down spins, 
each Hamiltonian is characterized by a Chern integer. 
While the sum of the Chern integers is zero due to time 
reversal symmetry, its difference defines a quantized spin 
Hall conduct ivity^'^^. The Z2 invariant is then given by 
half the difference of the Chern integers {i.e., the spin 
Chern number) modulo 2. Here we will choose another 
way and calculate the Z2 invariant directly for the KM 
model. We follow Fu and Kane^^ and briefly adapt the 
calculation of the invariant for the sake of completeness. 

To compute the Z2 invariant it is important to keep the 



full tight binding model; in this section, we use the no- 
tations of Ref. 11. We express the matrix T-Lk of Eq. (14) 
in terms of F matrices, 

5 

^fc = 5I'^-(fe)r" (49) 

a=l 

where the five F" matrices are given by F^ = (g) /, 
T"^ ^ tV ® /, F3 = (g cr^, F^ = (g and F^ = 
(g) (T^. The coefficients di and c?2 are essentially real 
and imaginary part of g (more precise, real and imaginary 
part of exp (— ifc^s) (7), respectively, and ds = 7. ^3 and 
(i4 are both zero. In addition to the five F° matrices, 
there are also their ten commutators F"'' = [F", F'']/(2i). 
Their anticommutators F^F'' + F^F" = 25ab obey the 
Clifford algebra. The parity operator is defined as 

V = t'' ®I = V^ . (50) 

Obviously, inversion V interchanges the sublattices (t) 
but not the spin (a). The time- reversal operator T is 
defined by 

r = i{I®ay)K, (51) 

where K denotes complex conjugation. The Dirac matri- 
ces are chosen to be even under VT, VTV^iVT)^^ — T°- , 
while the commutators are odd under VT. Note that all 
F° are odd under V and T except F^ which is even. Time- 
reversal and inversion symmetry imply that the product 
VT commutes with the Hamiltonian. The only time- 
reversal invariant points of the BZ which have to fulfill 
for a reciprocal lattice vector —Ti = Ti + G are given by 

r. = ^(^1^1+^262) (52) 

with rii = 0, 1. We define Fi as the F point of the 
BZ {i.e., k = (0,0)), r2 = ibi, T3 = i62, and 
r4 = -^{bi + 62) where the bj are the reciprocal lattice 
vectors of the honeycomb lattice as shown in Fig. 8. The 
latter three points are usually refered to as Af points. 
The Z2 invariants characterizing the occupied band are 
determined ^""^ by 

5, = -sign(di(r,)) . (53) 

The Z2 invariant = 0, 1 which distinguishes a topologi- 
cal band insulator in two dimensions from a conventional 
band insulator is given by the product of all Si, 

4 

(-ir = n'^- (54) 

i=l 

Since we can write di{k) = t{l + cos (fcoi) + cos (fca2)) 
and use aibj — 2n5ij, we find diiTi) = 3t, di{T2) — 
diiTs) = t, and rfi(r4) = -t. Thus Si ^ S2 ^ S3 = -1 
and Si = +1 which implies by virtue of Eq. (54) 

1/ = 1 . (55) 
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In fact, a non-zero Z2 invariant implies a topological band 
insulator phase provided there is an energy gap through- 
out the BZ. In particular, the argument with invariants 
is a single-particle picture argument^®. 

As the mean-field interacting Hamiltonian has been 
reduced to a single particle Hamiltonian (which leaves 
di unchanged) the argument based on Z2 invariants is 
applicable as long as the spin orbit gap is present. The 
band Hartree Fock approach has precisely shown that the 
Hubbard term behaves as a chemical potential which lies 
between the bands as long as \U\ < 2t. Thus we have 
substantiated our earlier statement that the TBI phase 
and the presence of edge modes will be stable up to a 
region beyond the perturbative regime. 




0.2 0.4 0.6 0.8 1 



A 

FIG. 9: (Color online) Numerical solution of the mean field 
equation (68). The behavior of Qx{^) is shown (t = 1). 



V. MOTT TRANSITION 

Now we will use the slave rotor mean field theory^^'^^ 
to find the transition where the charge degrees of freedom 



Here we introduced the renormalized KM spectrum for 
the spinon sector. 



form a Mott insulating state (and not a band insulator). = y {Qf I9I) + (Q/ 7) j (63) 



A. Self-Consistency equations 



and defined 



In contrast to the approximation in Appendix B we 
will use a more sophisticated approach ^5,76,80 order to 
find the transition to the Mott phase. We replace the 
phase field representing the 0(2) degree of freedom by a 
complex bosonic field X = e*^ which is constrained by 



\X{Tt 



1 



(56) 



The associated Lagrange multiplier is called p. The 
derivation of the Lagrangian and the decomposition of 
hopping and spin orbit terms is shown in Appendix C. 
The mean field parameters associated with the decom- 
position are given by 



Qf = (exp(-i%)) , 
for the hopping term and 



Q'f = (exp (-*%)) , 



(57) 



(58) 



(59) 



(60) 



for the spin-orbit term. Finally we find the imaginary 
time Green's function for the fields. 



Gfi - 



1 



and for the X fields, 
Gx 



(61) 



(62) 



ik^-Qx\gm + Q'x>^92{k) 



(64) 



with g2 from Eq. (19). Note that we consider only the 
half filled case here which allowed us to set fi — h = 
0. In the absence of spin orbit coupling, A — >■ 0, we 
find T^k ^ Qf \g\ and £,k — t- —Qx \g\ and recover, hence, 
the Green's function from Florens and Georges^^ (when 
setting e = — |g|). From there, we find directly the five 
self-consistency equations. 



U 



y I 

k jA2+4;7(e-min(a)) 



(65) 



In Eq. (65) we performed the evaluation of the Matsub- 
ara sum at zero temperature (see Appendix D) and we 
introduced the insulating gap 



min (^fc)) 



(66) 



While Ag is non-zero in the insulating phase, directly 
at the phase transition it will vanish since the rotors 
condense. Before we can use Eq. (65) to find the tran- 
sition line, we have to know the explicit form of (,k and 
hence of Qx and Q'x- The latter two mean field param- 
eters are determined by use of the (second and third) 
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self-consistency equations. We start with Qx- 



j = l a 



ha 



^ 5 +«;/?+) (/LVL) 



(67) 



-E 



Due to the lattice symmetry, the sum over the three near- 
est neighbors, j=i ^ j'^^* appears as a factor 3 in the final 
expression. Thus we find the mean field parameter 



\9? 



We have plotted Qx as a function of A in Fig. 9. 
In a similar way we proceed in order to find Q'x- 



(68) 



j acr' 



^5]7(i«-p-Kp-i/3-p+i/3+n 



X 



0.12 
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FIG. 10: (Color online) Numerical solution of —XQ'x (red 
curve) and —Q'x (blue curve) as a function of A (t = 1) is 
shown. 



which defines the transition line between TBI and the 
EMI phase as shown in the phase diagram Fig. 1. The 
sum over k' means that formally the lowest bound cor- 
responds to fc femin + V ^ Ij and femin IS associated 
with the minimum of ^fc. Hence, no divergence appears 
in the sum. A formal justification to cut the sum can be 
given by switching to "energy space" and considering the 
density of states. The same line of argument applies to 
Eqs. (73) and (75). 

As a last point we have to consider Q f and Q'^ and its 
behavior along the line C/c(A). Applying the same line of 
reasoning as for Qx we directly find 



-y 



27^ 



(69) 



Again the lattice symmetry is responsible for the fact that 



the sum over the next-nearest neighbors, 



can be 

replaced by a factor 6. Then the self-consistency equation 
reads 



Qa- =^E ^^^3<<y'ftafjaj 



(70) 



-|Qx(A)| 



We have plotted Q'x as a function of A in Fig. 10. With 
the knowledge of Qx and Q'x finally the rotor spectrum 
^fc of Eq. (64) is well defined and we can proceed with 
Eq. (65). If one moves towards the transition from the 
Mott insulator to the superfiuid phase of the rotors, the 
rotor gap Ag must close. It yields 



f/c(A) 



1 



H -2 



— E- 

V - min(a) 



(71) 



Qf = (xix, 



I] nn. 



Na 



Ly-MlyGx(fc,*i^„) 



1 y-|g| u 

Na ^ Qt 2y/U{p + ^k) ■ 



(72) 



Here denotes one of the three nearest-neighbor vec- 
tors. Along the transition line we have = and ob- 
tain 



)^fA^^ ^^ V I.9I 

'^^ ' etiVA if Ve'^-min(a) 



(73) 



It turns out that Q^- is a slowly varying function of A. 
We have plotted it in Fig. (11). 
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FIG. 11: (Color online) Numerical solution of QJ(A), i.e., Qj 
along the line Uc{X), as a function of A (t = 1) is shown. 
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FIG. 12: (Color online) Numerical solution of Qj{X) as a 
function oi \ {t = 1) is shown. 



The last self-consistency equation determines Q'r 



Q'f = {x:x, 



k 



u 

2NI 



ikS' 



(74) 



with S'^ being one of the six next nearest neighbor vec- 
tors. Thus we find Q'y along the Mott transition, 



P(A) = 



E 



(75) 



which we have plotted in Fig. 12. From Figs. 11 and 12 we 
see that Q'j: and behave similarly but is (roughly) 
three times larger than Q'jr. 





FIG. 13: (Color online) Left: Spectrum ±e{k^) of the KM 
model on a stripe geometry as explained in the text for A = 0.2 
and t — 1. Right: Spectrum ±E(fci,) of the spinon sec- 
tor for the same parameters at the critical line C/c(A = 0.2). 
The hopping and spin orbit amplitudes are renormalized with 
Q} = 0.68 and Qf = 0.26. The tildes refer to the finite stripe 
geometry in contrast to the spectra of the infinite system. 



B. Discussion 

In Fig. 13, we show the spectrum of the non-interacting 
KM model for A = 0.2 as well as the renormalized spinon 
spectrum for Uc{X — 0.2). The geometry we considered is 
a stripe with 14 unit cells in y-direction while the stripe 
is infinitely long in a;-direction. Here we see that interac- 
tions contribute through Qf and such that the bulk 
spin gap is decreased compared to the TBI phase. 

At the mean-field level, spin-charge separation will oc- 
cur and, while the charge is frozen in the Mott insulating 
state, the spin degrees of freedom will exhibit an Hamil- 
tonian reminiscent of the KM model. In particular, this 
implies the existence of gapless edge spinous. In this 
sense this gives rise to the Fractionalized TI mentioned 
in the introduction^^ or a "topological Mott insulator" 
as Pesin and Balents"*^ did. Note that the "topological 
Mott insulator" phase has a different meaning than in 
the work of Raghu et al.'^^ where the topological band 
insulator phase was caused by strong interactions. On 
the other hand, as already mentioned in the introduction, 
U(l) gauge fields associated with the lattice theory*^, see 
also Appendix E, cannot be ignored in two dimensions. 
In particular, to stabilize the Fractionalized TI beyond 
the mean-field level, one requires extra layers support- 
ing gapless spinous allowing to screen the gauge field^^ . 
(The stability of spinon excitations at the edges results 
from the fact that single-particle tunneling is suppressed 
as a result of Mott physics^^). This also implies that in 
the context of an isolated (single) honeycomb layer the 
Fractionalized TI is unstable to instanton proliferation 
and to easy plane Neel ordering^^'^'^. 

Let us assume that conditions are realized such that 
the Fractionalized TI is stable against gauge fluctuations. 
Then, we can be even more precise when focusing on the 
spinon bulk sector. We can write for the corresponding 
ground state wave function. 



1*/)= n f> 



0] 



(76) 
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i.e., the lower band — Sfc (see Eq. (63)) is completely 
filled. The explicit knowledge of allows us to calcu- 
late the expectation value of the ^-component of spin, 

k 

In the last line we have used Eq. (38). In the same way, we 
can easily check that (Sf) = 0. At the mean-field level, 
we can also calculate the spin-spin correlation functions 
{S^SJ) for i and j on the same sublattice or on different 
sublattices. We show the former case explicitly: 

{S+Sr) = //^VjVI*/) (78) 




where we can evaluate the last line numerically and find 
that the correlations decay to zero on very short dis- 
tances. In fact, when the distance \Ri — Rj\ reaches 
roughly four unit cells, the correlations are already 
smaller than 10~^. A similar calculation for both i and 
j on sublattice B as well as i and j on different sublat- 
tices reveals comparable results. We expect that the one- 
dimensional character of spinon excitations at the edges 
gives rise to power-law spin correlation fimctions. Other 
exotic spin liquid phases may be found in the vicinity of 
a Mott statc'^^. 

The only "hidden" order which seems to be present is 

reflected in (^E"/! X^fc 7 '^'I'^^'^^^fe I*/) being a rem- 
nant of the TBI phase; here we have introduced the cor- 

responding vector <^>l ^ (/^V ' /fet^ /fc/) : 

In the phase diagram of Fig. 1 and 2 the region above 
the TBI phase has to be handled with care for < A < 
0.1 i and is beyond the scope of this paper. In particular, 
in the absence of spin orbit coupling (A = 0), mean- 
field slave rotor techniques predict a Mott insulator with 
gapless spin excitations^^. In the limit A — >■ we recover 
the earlier result of Lee and Lee*", i.e., = 1.68t. On 
the other hand, in a recent QMC study it was shown*^ 
that the intermediate phase for A = is an RVB spin 
liquid in contrast to Refs. 80,82. 

While we have clarified the question what the effect of 
a Hubbard onsite interaction is one could also consider 



nearest and next nearest neighbor repulsion with ampli- 
tudes V\ and V2. Such a model in absence of spin orbit 
coupling was investigated by Raghu et al}^ . From the 
band Hartree Fock approach presented in Sec. IV we see, 
however, that the effect for small Vi and V2 is negligible. 
This is in correspondence with Ref. 43 where strong near- 
est and next nearest neighbor interactions are required to 
reach QSH phases while weak interactions leave the semi 
metal unchanged. As the intrinsic spin orbit interaction 
already opened a gap, V\ and V2 should be of the same 
order to affect the phase diagram. As large values of 
and V2 are in the model under consideration somewhat 
unphysical, we can conclude that additional weak nearest 
and next nearest neighbor interactions are negligible for 
the KM model. 



VI. CONCLUSION 

We have investigated the Kane-Mele model in the pres- 
ence of a Hubbard interaction as a paradigm for two- 
dimensional topological insulators with interactions. 

Using a mean-field procedure and arguments from the 
slave-rotor theory, we have shown that the TBI phase 
characterized by a Z2 topological invariant is stable 
(against Mott phases) up to moderate interactions which 
are beyond the perturbativc regime. The topological 
band insulator phase is separated from a Mott insulating 
region through ?7c(A). 

At the mean-field level, charge constituents become 
frozen in the Mott state while the spin constituents form 
a quantum spin liquid with gapless edge spinous (pre- 
serving the time-reversal symmetry). Even though this 
Fractionalized TI phase is unstable against gauge fluc- 
tuations in the isolated honeycomb lattice system^^'''^, 
the vicinity of other screening layers exhibiting gapless 
spinon excitations^^ would allow to stabilize such a phase 
of matter in (quasi-) two dimensional systems. For very 
large onsite interactions, the Fractionalized TI phase in- 
evitably turns into a SDW phase with XY ordering. 

For very weak spin orbit interactions, other insulating 
phases reminiscent of the "gapless Mott insulator" phase 
of Pesin and Balents^^ might exist. It remains an open 
question if such a possible phase might be connected with 
the expected spin liquid phase for A — )■ 0. 
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Appendix A: Band beisis 



Appendix B: Simple slave-rotor mean-field theory 



In this Appendix, wc present some supplementary ma- 
terial concerning the band basis which was introduced 
in Sec. IV. First, let us check. 



and t/T; 



that T^Tf = TfT^ = 1 



T^T^ = 1. We find the diagonal elements 

lo^ibl^ + |/3±P = 1 which can be easily checked. The off- 
diagonal elements are given by 



a*_a+ +P*_P+= J\f-Af+{d*_d+ + l) 



. 



(Al) 



We will further need the following expressions, 



AA+A/l = 



2s 



and 



2(£2 ± js) 



(A2) 



In the above transformations the limit g ^ should be 
handled with care as the original eigenvectors diverge. 
This is a consequence of the fact that the matrix Hk is 
already diagonal for t = 0. To consider the case with 
A = 0, we recover the transformation matrix (9) from 
Section II, 



lim Tt = hm T, = Tq . 

7^-0 7^-0 



(A3) 



In order to prepare the following part of the Sec- 
tion, we have to show explicitly, that X]fcl'^±(^)P = 
Efe l^±(fe)P = One can show that |a±|2 = |/3^|2 

(which implies |a+p + |a_p = 1) and it is sufficient 
to consider \a±\'^ in the following. First we show that 
|a±(A;)p + |a±(— — 1 and then we argue that 
the result follows directly. By using |(7(— fc)| — |5(fe)|, 
e(— fc) = e(fc), and j{—k) = — 7(A;) we could calculate 
|«±(fc)P + |q;±(— = 1 exphcitly. This is not neces- 
sary since the only thing we have to show is 



\a.{k)f = \a+{-k)\' 



(A4) 



By looking at the definition of a± and using 7(— fc) = 
— 7(fe), Eq. (A4) turns out to be correct. Now we have to 
divide the BZ into two parts, e.g. as follows: BZ=/Ci U/C2 
with /Ci = {fc e [-^,0] X [-^,0] and fc € [0, ^] x 

[— ^^,0[} while /C2 contains the remainder of the BZ 
such that /Ci n /C2 = 0. This ensures that fe € /Ci implies 
—k € IC2 and vice versa. 
Hence we can split the sum over the BZ as 

^ |a±(fe)P = ^ |a±(fc)|^+ l"±WI' 

feeBZ kelCi 

= (|a±(fe)|2 + |a±(-fe)n 
keiCi 

= ^ 1 = 7Va/2 . (A5) 
fee/Ci 



In this Appendix, we perform a simple mean field de- 
composition for the KM model with Hubbard interac- 
tions which is rewritten in spinous and rotors. Starting 
from Eqs. (46) and (47), adding the rewritten interaction 
term Hj = - Lf , we assume that state vectors in 

the Hilbert space should have the form 1^*) = l'^ f) l^'e)- 
Decoupling the rotor and fermion variables and treating 
the constraint by introducing a Lagrange multiplier h, we 
obtain effective Hamiltonians for spinon and rotor sector: 



{ij)iT <Siij'> og' 

-(M+/i)E/Ua> (Bi) 



- E cos (^'^) + E ^ cos 

{ij> «U» 



+E 



hL, 



(B2) 



In the spinon sector we recover the original KM model 
expressed in f„ operators rather than c„ operators while 
the rotor sector corresponds to a quantum- XY like model 
for the phase variables. The effective amplitudes are de- 
termined by the following self-consistent equations. 



t|f = t(cos(%)), 
= A(cos(^,,)). 



-T-eff 



(B3) 

(B4) 

AE(^^'i^-'/^/.>')/ • (B6) 



The expectation values are taken with respect to \^ or 
Ivfe), respectively. In the following, we show the main 
steps in finding the self consistency equations and the 
phase transition within this simple approach. In fact, 
it is the analogous calculation for the honeycomb lattice 
to Sec. HI. B. 2. and 3. of Ref. 76 for cubic lattices. We 
restrict ourselves to the case A = 0; the Hamiltonian "H-^ 
reads in momentum space 



T-i^ = - E ( fka ! /fecT ) 

T^^ ^\ Zg* 




(B7) 

Here we decomposed the expectation value simply as 
(cos {9i — 9j)) PS (cos = Z. By diagonalization of the 
matrix in Eq. (B7) we obtain the upper and lower band, 



e± = ±Z\g\ - {11 + h) . 



(B8) 
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Wc stress again that our constraint is differently chosen 
compared to Ref. 76 where it is given by fjafia ~L = 
1. In the constraint used in this paper the L term appears 
with a positive sign. The difference is due to the different 
definition of the fcrmions in term of spinous and rotors, 
see Eq. (6) of Ref. 76. Now, the Lagrange multiplier h is 
defined from the following constraint equation: 



(L), = -^((/U.)/-0 • 



(B9) 



To proceed further, we treat "H^, which corresponds to a 
quantum XY model, at the mean-field level. The applied 
approximation, cos {Oij) « 2 cos (6'j)(cos (0^)) — const., re- 
duces the rotor Hamiltonian to a mean-field Hamiltonian 
of independent sites: 



(BIO) 



Here the coupling constant K is given by 



a t /-b 



(Bll) 



As long as we are in the rotor-condensed phase, we can 
assume (cos^j) = (cos 6) which allows us to evaluate K, 



K = 2(cos^)5^-^(-i)5^(5A)(/LVL ) 

= 4(cos^)^^-5|^(/LVL > (B12) 
= 2^^-|3|(cos0)s2e-(l/2)(cose) . 



The matrix elements in the previous equation are easily 
calculated using the definitions of Sec. IV, 



{fkJfka ) — ' 



-al/3-(4t^fet) 



_ 9_ 

2 \a\ ' 

~ 2|g| • 



(B13) 



We define the half-bandwith D = 3t for the honeycomb 
lattice, and find numerically the result 



^E</"c.^4 ) = l^"(l/2)|^1.57 



(B14) 



Now let us consider Eq. (B9), 

(i) = -^((/U.)-0 



= 0, 



(B15) 

where we assumed that the lower band is completely filled 
while the upper band is empty. Also, since we still assume 
the half filled case, Eq. (40) of Ref. 76 is easily 

r^=^E(/t/-) = ^E^ = ^ (B16) 

This is equivalent to set fj, + h = 0. Thus, we can also 
introduce /Xo(n) which is defined as: 

h + fJ, half filling 



A*o("-) 







(B17) 



Similar to Ref. 76 we obtain the Green's function 

Gfl{k,iuJn)~^ = ILOn - ZSk ■ (B18) 

As a last step {cos0) has to be calculated. Following 
Ref. 76 we calculate it in first order perturbation theory 
in K. We start with (cos^) = (V'l^'l cos^|V';^^) where 

{I\c0s9\ln) 



W3 = iw + 



\i) . 



(B19) 



To first order in K only the "mixed" element contributes: 

(/|cos0|UP 



(cos 61) = 2KY^ 



E, 



(B20) 



The energies arc given by Ei = l/{2U)(UL + h) -|-const. 
and the matrix elements by ||(^|e*^ + e~*^ |Z„) p = 
j(<5i,z„_i + <^;,i„+i)- Altogether we find the result 

2K 



(cose) = 



u 



(B21) 



which is in agreement with Ref. 76 at half filling. Then 
we substitute {cos 9) in Eq. (B12) and obtain finally 

?7e = -4£(l/2)=4|£(l/2)| . (B22) 

By means of Eq. (B14) we find the phase transition at 

J7~~6.30i, (B23) 

which should be considered as the correct result in 
d = oo dimensions. In this Appendix, we used a sim- 
ple mean-field approximation with the severe restriction 

Z = {cos (dij)). This approximation might be justified 
in large dimensions. Therefore, we can assume that the 
result U?° is exact in = oo. 
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Appendix C: Derivation of Green's functions 

In this Appendix, we pedagogically show all the rel- 
evant steps starting from the slave rotor Hamiltonian 
to the Green's functions. We have omitted this part in 
Sec. V for the sake of clarity. The slave rotor Hamiltonian 
reads 

^ = -tEE(/"-^-/i- e-^^«+h.c.) 

(ij) <^ 




where we still have to fulfill the constraint Eq. (45) with 



the Lagrange multiplier h. Then, the action is built from 

Sq= dT[-iLdr9 + U + f'^drf] , (C2) 
Jo 

where the first two terms correspond to the Legendre 
transform between H and jC and we are switching from 

phase and angular momentum operator (9, L) to fields 
{0,drO). Here L and drO are related as follows. 



= -wr , (C3) 



which yields L = {i/U) drO. We obtain the action 



^0 =iy^ E (^x - M + /li) /ia + ^ E i^-^i + ' 



E 



2U 



{ij) \ a- 



^ E (E^^*^<-'/-/^-e"^'"^) 



(C4) 



Now we have the choice to decompose the hopping and 

spin orbit terms either in a standard way (as Florens and 
Georges did*"^) or in a more elaborate way (as Lee and Lee 
did^*^) to obtain an effective theory. Since we are mainly 
interested in the transition line to the Mott phase we 
restrict ourselves to the first way for the moment. Thus 
we will use again the decomposition a/S « {a)l3 + a{j3) — 
{a) {13} with 

E./r/j-., {aij) =Qx, 
Pij= exp{-ieij), {I3ij) =Qf, 



for the hopping term and 



= cxp(-i6'y), 



K) = Q'x , 



for the spin-orbit term. Then we replace the exponen- 
tials exp(i^j) by complex bosonic fields X{t) which are 
constraint via jX^p = 1. This constraint is imposed by 
a complex Lagrange multiplier p^. Then the decomposed 
Lagrangian has the form 

S = S' + S" + S'" , (C5) 

where S' contains the hopping term, S" the spin orbit 



term, and S'" the other terms. S' is given by 



Jo 



dT 



\t\QxY.^tXj +C.C. 



\AQf E E ^tafj'^ + c.c. + \t\ E QiQx] (C6) 



Jo 



dT 



jCx+jCf + ... 



The second term S" is given by 



S" : 



dT 



XQ'x E ^tXi + c.c. 

{{ij)) 



+ ^Q'f E E*^*^<-'/i'-/^-'+^ E Wx 

{{tj)) crcr' ■^ij:S> 



: / dT[£.'x+L'f 

Jo 



+ ... . 



(C7) 



Introducing the X-fields changes the Hubbard term such 

as 

{dr9i + hif = [{idr + hi) X*] [{-idr + hi) Xi] , (C8) 
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and the term S'" becomes the form 



2U 

+ ^ + ^ (a, - M + /li) + . 

^0 

(C9) 

Here and in the previous equations the eUipsis corre- 
sponds to the other terms which are independent of f„ 
and X. The Fourier-transform of Cx and Cf leads in the 
(u, I) basis to the bands obtained earher: 

jCx = QxY,{-9{k))Xl*Xl + c.c. 

k 

= QxJ2(-\9\)Xl''xl + \g\XrX]i , (CIO) 
fe 

= Q/^(-5(fe))/feV/L+c.c. 

ka 

= QfJ2(-\9\)fL*fL + \9\fVfL-{cn) 

her 

For the X-part of the spin orbit term C'-^ we find the 
following expression: 



c'x = Q'x E ^92{k) (xrx£ + X^X^ 



Q'xY.^92{k) (xi^xi+xrxi) , 



(C12) 



where (72 (fc) is the; usual next-nearest neighbor hopping 
contribution as defined in Eq. (19). The last term which 
must be transformed into momentum space is C'^ which 
clearly produces the 7-term. Therefore we will add >C/ to 
£j in the {f^,f^) basis and then transform both terms 
to the (/^, /^) basis as we did with the original bands of 
the KM model: 

Cf + 4 = E Qf {-9{k)fVfL - gikffVfL) + 

ka 

^aa'Q'fl {fkcr^fkcr' ~ fka /fea') 

kaa' 

= X/ f'^<^ ff^<^ fkcT*fkcr ' 



ka 



(C13) 



Here we have introduced the renormalized KM spectrum 
for the spinon sector, 



^fc = \/{Qf \9\y + {Q'fiY ■ (CM) 



Finally we find the imaginary time Green's function for 
the fl fields, 



G 



1 



fi 



iuJn — Sfc ' 



and for the X fields. 



G 



1 



X 



(C15) 



(C16) 



where we defined 

Cfe = -Qx\9{k)\ + Q'xX92{k) . (C17) 

Appendix D: Matsubara Sum 

In the self-consistency equations (65), (72), and (74) 
we had to evaluate the following Matsubara sum: 



^^Gx{k,ii^n) = ?E 



1 



(Dl) 



1 



u 

^ {iv^ + A){-Wn + A) 

where A = \/U{p + ^k)- By taking the corresponding 
contour, we find 



= 



c {z + A){z-A) 

= E 



dz 



(D2) 



-2-Ki 



nB{A) nB{-A) 



2A 



2A 



where nB{z) = (exp (/3z) — 1) ^ is the Bose function. The 
last equation then implies: 



^^2ii^n + A){-iUr, + A) 

^ coth {13 A/2) T^o ^ 
2A ^ 2A ' 

At zero temperature we find the result 

U 



nejA) -nej-A) 
2A 
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Appendix E: Dynamical Gauge Field 

In this Appendix, we decompose the action Eq. (C4) 
in a different way compared to Appendix C and show, 
how the effective field theory and the gauge field Uij 
emerge. Following Lee and Lee^° we decompose the 
hopping term via Hubbard Stratonovich decomposition, 
«i.=E./ri^.andA,=e»''^ 

j dr]ijdr]*jdr]jidr]*i x 

■^g-i[\Vi3\^ + \V3'\^-Vij/3i3-VijO!ij-v*il3ji-Vjiaji] ^g]^^ 
~ P 

The same equation holds for a' and ^' in order to decou- 
ple the spin orbit term. 



In Eq. (El), £ is given by At times the hopping or the 
spin orbit amplitude, respectively. We further follow 
Lee and Lee and change the variables of integration by 
V^J = |x,;,k-'"-+^(4---) andr?,, = |x,,k-™-+'(4--). 
Note that rjij and r]ji are independent complex variables, 
and hence Wij and a^- are independent and necessary. 
At this point, we replace again exp(i0j) by the bosonic 
Xi-field with the constraint \Xi\'^ = 1 imposed by the La- 
grange multiplier pi . Then we find the action which coin- 
cides with Eq. (4) of Ref. 80 apart from two terms coming 
from the spin orbit interaction and the slightly different 
notation of the rotor variables. We replace the variables 
by their saddle point values plus fluctuations (see for de- 
tails Ref. 80) , neglect the massive modes and we can in- 
tegrate out the Pi field to restore the 6 field. Finally we 
obtain the effective Lagrangian (similar to Ref. 80) : 



P'ij = exp(-i%). 



io- i {ij),cr {ij) 

(E2) 



Here aj" and are the temporal and spatial gauge fields 
coming from the fluctuations from hi and Uij, respec- 
tively, hi is the Lagrange multiplier associated with the 
global constraint introduced earlier. The quantities with 
tildes are the saddle point values and are identical to the 
mean-field parameters which we have evaluated in Sec. V. 
We notice that both spinous and rotors couple to the U(l) 
gauge field. Since we assume weak gauge fiuctuations, we 
take the saddle point approximation, i.e., = 0; we re- 
cover for the spinous the same terms which resulted in 
Sec. V in the renormalized KM spectrum (xjj — >■ Qx and 



Q'x)- We should also mention that the spinous 
are still coupled to the gauge field through the first term 
in Eq. (E2) which contains ^ f*^aj fi„. In principle, al- 
though we have set aj =0, we could allow for small de- 
viations and expand exp(ia^) « l + ia-jt; thus the spinous 
couple to both temporal and spatial gauge fields. On the 
other hand, we know that the rotors are gapped in the 
Mott phase and can be integrated out. This generates 
the Maxwellian term (see e.g. Refs. 72, 73, 80). 
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